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I . Introduction 

This is the final report of researcli conducted since 
1976 under NASA Grant NSG-1243 under the general title "Rain 
Scavenging of Solid Rocket Exliaust Clouds". The requirement 
addressed by the work is that of assessing the environmental 
impact of Solid Rocket Motor (SRM) exl\aust products dis- 
charged into the free air stream upon the launcliing of space 
vehicles tl\at depend upon SRM boosters to obtain large thrust. 
These include the Titan series and especially tlae Space Shuttle. 

The exhaust product of greatest concern is HCl gas, 
of which the Space Shuttle boosters generate and discharge 
128,900 kg to the air below 10 km (troposphere). In addition 
142,543 kg of H^O and 174,900 kg of AI 2 O 2 are emitted to the 
troposphere from the SRM s in each Shuttle launch. The AljO^ 
appears as particles suitable for heterogeneous nucleation of 
HCl^^ (hydrochloric acid) which under frequently occurring 
atmospheric conditions may form a highly acidic rain capable 
of dajnaging property and crops and of impacting upon the 
health of human and animal populations. 

The meteorological assessment of this problem has 
numerous aspects. The present work has addressed two of these, 
najnely, (a) the cloud processes leading to the formation of acid 
rain and the concentration of the acid that then reaches the 
ground, and (b) the atmospheric situations that lead to the 
production of cloud and rain at and near' a launch site, and the 
prediction of weather conditions that may permit or prohibit 
a launch operation. 


1. 


In Section II of this report we present our analysis of 
the heterogeneous condensation/evaporation of HCl and H 2 O under 
conditions found in Titan III exl\aust clouds (’’ground cloud”) 
some 90 sec after launcli at about 1 km altitude. TI\is provides 
basic information that should be used in a cloud/rain micro- 
piiysical model to predict rainfall occurrence and acid 
concentration. 

Section III presents a numerical method for use in 
generating weather predictions by moans of our 3-D Mesoscale 
Model (Hsu, 1979) . 

II. The CQ-CQndetAsation/Evapoi‘ation of HCl and H,0 

«• 

The free energy change, AF, of a system is a measure of 
the tendency of that system to progress from one thermodynamic 
f vato to another. For the case of condensation/evaporation of i 
different vapors on wettable particles, the general expression is 

AF * ^cf’Ca" - r “) - ^ (h. kT In S>’) (1) 

p i 1 1 

whore 

o' is the surface enorg)'' per unit area of the droplet 
surface 

a is the radius of the droplet 

r^^ is the radius of the nucleating particle 

is the number of i-molecules condensed on the drop 

k is Boltzmann’s constant 

T is the absoliite temperature of the droplet surface 

$.* is the saturation ratio of vapor i witli respect to a 

^ flat surface of bulk solution at the same concentration 
as the droplet 


w 


The present purpose is to explore the application of 
this basic thermodynamic statement to the case of the co-condensa- 
tion/evaporation of H^O and HCl vapors on wettable particles 
in the open air. For this case, let i » 1 specify 11,0 and i * 2 
HCl, aJid lot the HCl molofraction, x,, express the solution con- 
centration. By definition, then, the molality, M ® 55. 5 
Xo/Cl - x^l =55.5 f(x^). 

The drop radius is 


- n.m, •: 1 / - 

“ = [47 ^ 1'^ 


whore 


energy, 


m, is the mass of a molecule of H^O 
p’ is the solution density 
Hq is Avogadro’ s niimber 

B is M-,/M, and M. is the nralecular weight of i. 

*■ X X 


Empirical expressions are used for the solution surface 


o' = 75.728 - Q.1535CT-273.16) - 10.575fCx„) 


and the solution density 


p' « 1.0 - 0 . 72158 fCx.,) 


The saturation ratio for ead\ vapor with respect to 
a flat surface of the solution is given by 


Si' « P./P3 ^^'^^Cx2,T) 
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where 

is the environmental partial pressure of i 
and 

P. Cxj.T) is the equilibria, vapor pressure of i 
over a flat surface of solution of concentration 
and temperature T. 

For T = 288° K, the values of S^' and S^' as functions of X 2 are shown 

-2 

in Pig, 1. Note that S,' > 1 only for X 2 < 8.85 x 10 , whereas S^’ > 1 

_o „ , , 

only for x., > 8.15 x 10 hence at 288 K, the vapors are both 

supersaturated with respect to the bulk solution only in the 

-2 -2 

narrow HCl molefraction range of 8.15 x 10 ^ 1 . 

By means of ( 2 ) , (3) , (4) and (5) , AF may now be expressed 
in terms of the six variables P., P_, T, n. , r , and x_. For any 

i £ 1 P ^ 

particular case, the environmental values of P^^, P., and T must be 
specified, thus AF may be expressed for such a case in terms of r , 

r 

n^ and x, . 

The total number of molecules, n^, required to form a 

monolayer of solution on a particle may be estimated as a function 

of r and x„ as follows. If the particle and its coating of solution 
P 2 

is spherical, then is given by 

r 2 

n = h, + n, = 4 (1 7 ^ ) C 6 &) 

s X - ^ 

s 

r 2 

also = 4(1 + 7 ^ W (1 + fCXj}) ( 6 ^) 

^s 

where r is the mean molecular radius of a ’’molecule" of solution, 
s 
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defined by 


OF FGvli 




[Cl - xj V, + X, V, ]} 


1/3 


where and v^ are respectively the molecular volumes of H,0 
HCl in solution. In general, the jnolecular volume of species 
be written 


and 
i may 


Vi = Mj^/p’no 

giving 

V, = 2.9914 X 10““^ C 
1 P ' 

and V 2 = ev^. 

• -S -1 

Values of r for the concentration range l.S x 10 < x., < 6.4306 x 10 

are given in Table 1. 


Table 1. 

Values of r for 
range 10"“^ < M < 

HCl in the 
,aq 

10“ 

concentration 


^2 

M 

^1 

10^-rgCm 

1.8 X 

10"^ 

10"^ 

0.999982 

1.9256 

1.8 X 

10" 

10-2 

0.00082 

1.9259 

l.S X 

10“^ 

10"^ 

0.9982 

1.9269 

1.77 X 

10 “ 

10° 

0.9823 

1.9373 

1.5266 

X 10"^ 

10^ 

0.84734 

2.0212 

6.4306 

X 10“^ 

10'2 

0.35694 

2.2794 
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A reasonable value for r^ in the concentration range of 
8 X 10**^ £ X 2 £9 X 10”^ is 2 X lO"® aiu Using this value in 
( 6 b), nj^ may be calculated, ajvi. AF may be determined as a 
function of r^ and x, for specified environmental conditions. 

Taking environmental conditions as found in the stabilised 

ground cloud generated by a Titan III launch: T = 298.16“K, 

0 2 

= 23582 dynes/cm“, = 104.5 dynes/cm at time t *= 90 sec., 
the map of AF in (X«,r ) coordinates (Fig. 2) is constructed. 

£. p 

The contours of the AF surface in r , x_ coordinates show 

p i 

a definite "saddle point", t, which for the specified environ- 

-6 -2 

mental conditions, occurs at r = 1.3 x 10 cm, x» = 8.8 x 10 . 

p i 

This point is analogous to the critical point defined for single 
vapor heterogeneous nucleation (see, e.g. , Byers, 1965, Chap. 2). 

The surface thus defined is equivalent to the fx-ee energy surface 
AG (n^,njj) for the H.,0 - H. 3 S 0 ^ system that has been discussed by 
Reiss (1950), Kiang and Stauffer (1973), Haraill (1975) and Hamill, 
et al (1977) . Reiss (1950) showed that, when the surface energy 
tem\ is included in the free energy expression, the free energy 
surface AG is saddle-shaped with a saddle point defined by 

^ Al 

(AG) = 0 and (AG) = 0. The saddle shape was also found to 
A B 

be present under stratospheric conditions for the H 2 O - ^ 280 ^ system 
by Hamill, et al (1977) . 

Several features of the AF (.x.,,r ) surface (Figure 2) merit 

im p 

-2 

discussion. At constant X 2 = 8,8x10 , {M - 5.355), the AF values rise 

*8 

gradually with increasing particle site in the range 10 ” cm j< r^ £ 

t T 1 

1.3 X 10 ” cm, reaching a maximum of 25 x lO” erg at the latter 
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Dry Pariicle Radius, rp, ia cm 




size (saddle point). AF then decreases to 0 at r == 2.67 x 10“^ cm 

P 

and decreases sharply to large negative values for larger sizes. 

_2 

At higher and lower concentrations, e.g., X 2 ^ 13 x 10 and 
-2 

X 2 4 X 10 , the AF surface rises sharply as r^ increases above 

“6 

10 " cm. Thus the growth region for HCl aq. droplets lies between 

-2 

nearly vertical ''canyon" walls at X 2 ~ 5.8 x 10 and X 2 ~ 11.7 x 10 

for nuclei of size r " 2 x 10~^cm and larger. The bottom of the 

P 

"canyon" is relatively broad and flat with the locus of minima lying 
-2 

near X 2 = 8.8 x 10 . 

The profile of the AF (x_,r ) surface taken at r = 10~^cm 

^ P P 

(Figure 3) shows the shape of the "canyon." In addition the 

values of the respective components of AF are shown for r^ = 10~^cm 

as a matter of interest. This diagram necessarily is discontinuous 

at AF = 0. The terms represented are: 

2 2 

AFc = 4na' (a - r ) - constant 

6 p 

AF^ = n^ k r in S^' 

AF 2 = A 2 k T In S 2 ' 

AF' = AF^ + AF 2 
AF = AF' + AF 

r 

AF 2 is affected by two factors: (a) as X 2 increases, ^2 also must 

increase, and (b) as X 2 increases S,' must decrease causing In S 2 ' 

to go from positive to negative values as S,' decreases through 1.0. 

-2 

The result of this is a minimum in AF 2 near X 2 = 4.5 x 10 , and a 

-2 

sign reversal near X 2 = 8.8 x 10 . No minimum is found for AF^^ 
because both n^ and decrease as X 2 increases. Inasmuch as 
droplet growth cannot proceed unless both vapors are saturated, the 
curves AF^ and AF 2 indicate that the region for droplet growth is in 
fact much narrower than the AF "canyon." This is true because, if 


o 



T n 




either vapor is undersaturated, a droplet in that vapor must yield 

the undersaturated species by evaporation because the mass diffusion 

is proportional to P^(S/-1). The approximate range for nucleation 

-2 -2 

to occur is 8.15 X 10 X 2 5 ^ 8.85 x 10 . The minimum value of 

-2 

AF occurs at X2 = 8.82 X 10 , thus this is the most probable HCl 

-5 

concentration for a particle size r^ = 10 cm under the specified 
conditions. 

It is clear that AF (X 2 , r^) varies with T, and P 2 . For 
the purposes of NASA, the maps of AF(x 2 »ip) for the various environ- 
mental conditions that may be encountered at different launch sites 
and in all seasons should be computed. Particularly the changes 
imposed upon the system by the increase of solid rocket booster 
capacity required for the Space Shuttle should be more completely 
evaluated. 
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An Explicit Mixed Numerical Method for 
Atmospheric Models 


by 


H-m Hsu 

Assistant Professor of Atmospheric Science 


University of Wisconsin, Milv^auKee 


Abstract 


An explicit, mixed numerical method has been developed for atmospheric 
models. In a set of physical equations, the forward finite-difference 
scheme is applied for the time tendency terras, upstream for the advection 
terms, and central for other terms. For either the shallow-water equations 
in one or tw dimensions or the primitive equations in three dimensions, 
the mixed method is conditionally stable and shows much better accuracy 
than tliat of the pure forward-upstream method. It is also shown that the 
traditional CFL condition is only a special case of the stability conditions 
revealed in this study. 
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1, Introduction 

Numerical simulation has become a more and more important method 
to reveal atmosphere processes, because the traditional analytic 
method has frequently failed to provide solutions from complex systems 
of partial differential equations which describe atmospheric phen- 
omena. Different numerical schemes are used to approximate such 
systems from differential form to difference form. Due to the limitation 
on computer resources^ economy and accuracy of the numerical scheme 
should be simultaneously considered. 

The central finite-difference schane has been widely applied in 
atmospheric numerical models, and recently became so popular that one 
tends only to empliasize its advantage in higher-order accuracy of 
solution and de’-emphasize its bad performance in actual calculations. 

In fact, this scheme not only generates erroneous small-scale pertur- 
bations which gradually distort the model results, but also provides ’ 
tw separate solutions for odd and even time-steps during the numerical 
integration. To avoid those errors . artificial space and time smoothers 
are necessary and have to be implemented into the computational algo- 
rithm. Among other explicit schemes, a lower-order scheme \d.thout 
any smoother may become an adequate alternative. After all, the 
central-differencing scheme (a three- time-level scheme) requires more 
programming efforts and greater computer x'esources than lo'tver-order 
ones (i.e. t^-xr- time-level schemes). 

The explicit t\70- time-level schemes used in atmospheric models have 
been described and summarized by Thompson (1) , Mesinger and Arakam (2) 
and Haltiner ’and Williams (3). For a single linear advection equation, 
the for\>xird-in-time and upstream-in-space scheme has been proved to 
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be stable, whereas the fonrard-in-time and. Qaatrai-inrsoace Is nnatable,. 
If a more complicated system than a simple advection equation is approx- 
imated by a mixed method, the stability of the approximate system cannot 
be safely determined without careful analysis. The mixed method 
considered here consists of the for^^;ard scheme for time tendency terms, 
the upstream scheme for advective terms, and the central scheme for 
other terms in the system. 

The purpose of this paper is to demonstrate that the mixed method 
is conditionally stable for both linearized systems of shallo\^^ater 
equations and primitive equations. 

2. Linear shallow-i^ter system 

Tlie shallo\^r\^ater system is the simplest primitive equation system 
for an incanpressible, hydrostatic, adiabatic and frictionless fluid. 
Kasahara (4) applied the central finite-difference scheme to a one- 
dimensional system with ti«) different staggered grid-nets, and analyzed 
the numerical stabilities. From his study he indicated that stability 
analysis should not be performed separately for every physical factor 
in the system, but for the entire system instead. Same shallow^^^ater 
equations will be adopted here to illustrate the numerical stability 
of the mixed method for those equations. 

The one-dimensional shallo\>rwater equations are 


3u 

3t 


, fi 9u , 3h _ 

+ U 0 


(l.a) 


IS. 


OF POOR QUAUtV, 


3h 

Tt 


u|^+ H 



(l.b) 


Symbols are defined in the Appendix. Basic properties of the mixed 
method may be revealed as follows, 
a. Stability 

Let denote the finite-difference approximation to u(t,x) S 
u(nAt,jlAx) and define h^. in a similar manner. The difference equations 
of (1) are 


n+1 


u. 


” u 


n 


n 


n 


.. 0 

Ax 


li 


g- 


n 

S .+1 


- 


‘a-1 _ 


= 0 


2Ax 


(2.a) 


and 


. n+1 . n 

hi h + u 


, n _ , n , n „ , n 

“ II - 0 


at Ax 2Ax (2.b) 

if U > 0. Sane set of finite difference equations is separately constructed 
by Bro\<m i?nd Pandolfo 5 j and stability is analyzed in an uncoupled approach. 


The solutions are assumed to have the following form 



Substitute (3) into (2), and an amplification matrix is obtained 
(Richtmyer and Morton, (6)). 
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G = 


~ i4) 


\l) - it|) , 


(4) 


•v 


where '1> = 1“U"^(1'" cos aAx), 

(5) 

({) = 0 sin aAx, 

(6) 

= S S 

(7) 

Ati 

and sin aAx. 

(8) 

The eigenvalues of the amplification matrix, G 


are X = - i (c{) 

(9) 

and X,.^ = - i (ct) . 

(10) 


While X is the physical mode of the solution, X^,. is the computational 
mode. The squared absolute value of eigenvalue is 
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and 


2 = 1 - 


1 - BAt + AAt" 


where A = 2(-^)^ (1-cos a Ax) + sin^ a Ax 


U . 2 

+ 2 ^ sin aAx, 


(Ax)' 


U 


B = 2(-^) (1-cos aAx), 




=/ir 


( 11 ) 

( 12 ) 


(13) 

(14) 


The von Neuman condition for stability (Richtmyer and Morton, (6)) 
requires that 


xr < 1. 


(15) 


Combining (11) and (15), a stability criterion is arrived, 

it <-|- ' (16) 

This concludes that the mixed numerical method is conditionally stable 
for a linearized one-dimensional shallown^ter system. 
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Theoretically, the stability bound emu be reached asymptotically, 


At,Y — 


U Ax 

(U + Cg)^ 


(17) 


as the wave number of a single Fourier component approaches to zero.'. 
Obviously, in term of velocity U, there is a maximum of the asymptotic 
stability bound for diffeirent depth of fluid. Figure 1 shows the 
^7)ariations of At,u for different depths of fluid and horizontal incre- 
ments under certain wind regime. If there is no surface gravity wave 
(H = 0) , the CFL (Courant-Friedrichs -Lewy) condition is met, and Atj;- 
decreases monotonically as U increases (Fig. la). For slow-moving 
waves. At,., increases with increasing velocity, then decreases after it 
passes its maximum (Fig. lb and c) . Finally, ALu will reach zero as 
0 approaches to infinite. For fast-moving \r;aves, the maximum of U 
is beyond 20 m/sec, and At,u decreases with decreasing velocity mono- 
tonically to zero (Fig. Id and le). 

The CFL Condition (Fig. la) is only a special case in the present 
results which give opposite conditions (e.g. Fig, Id and le) to the CFL 
conditions in certain situations (H = lOOra and lOOOra, respectively). 

It is also observed that for the shallower fluid At,u is dominated by the 
fluid velocity, and Hot the deeper fluid At^ primarily depends on the 
speed of surface gravity wave. 

The maxima for At,j^ and 0 can easily be determined, and they are 


and 


U = C 
max 


At = 


max 


g’ 

1 Ax 


(18) 

(19) 


Hence, it is concluded that Ac^^^ is inversely proportional to the 
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speed of surface gravity \<taves. In other words, the deeper the fluid, 
the shorter the tirae-<itep. 

Tliere are special cases and they may be described as follows, 
Case 1: If Cg = 0,> (11) becomes 


x|^ = 1-2 (^) (1 - COS oAx) At 


+ 2 (1 - COS aAx) At^. 


( 20 ) 


Thus 


At 


Ax 


( 21 ) 


is the condi tional stability criterion for the pure fonvard-upstream 
method applied to the simple advection equation. For simplicity this 
method is called pure method, which heavily criticised by Moelnkamp 
V) because of its highly dissipative character. 

Case 2: If U = 0, (11) becomes 


"" ^ ^ uAk)^. 


( 22 ) 


I I 2 * • • 

X is alwys greater than unity, and this method is absolutely unstable. 
Hence, U is not allo\H?ed to vanish, 
b . Damping factor 

The absolute value of the eigenvalue for the physical mode is a good 
indicator to ccxnpare how much the original wave amplitudes are reduced 
by the tiruncation errors due to different finite-difference methods. For 
the mixed method, jxj^ is described in (11). For the pure method, the 
third terms in both (2a) and (2b) are approximated by the upstream-in- 
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space scheme rather than by the central-in-space one, then we have 

Ixp = 1 - 2 [l -(U+Cg) (U+Cg) (-^) (1-COS aix), (23) 

following the same mathenatical procedure to obtain (11). The stability 
criterion is immediately provided, 

At (24) 

• U.Cg 

To shc«^ the comparison of the damping factors between these methods, ’ 
jxj^ is plotted against At for nine waves. In general, the time step 
for a stable calculation is more restrictive for the mixed method (Fig. 2) 
than that for the pure method (Fig. 3). Ho^^ver, |xj ^ gives remarkably 
less damping for the mixed method than that for the pure method. For 
example, the least X for 5Ax wave is greater than 99.5 5% for the mixed 
method, and is about 75% for the pure method. The high accuracy of the 
mixed method is obvious. Another interesting point is that the time step 
decreases mth increasing wavelength for the mixed method, ^dnile it is a 
constant for the pure method. The indication is that an accurate repre- 
sentation of requires short time step and the truncation errors due 
to fore^>;ard-in-tirae difference scheme is greatly reduced. Over all, it 
is concluded that the mixed method is* more accurate even though relatively 
small time increments are necessary. Also, the mixed method demands the 
same programing efforts and canputer resources as the pure method does, 
but much less than the leap-frog method, 
c. Phase speed 

Tlie false computational dispersion associated \vLth a finite-difference 
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method usually distorts the true solution of the problem. The acceleration 
and retardation of the approximate solution from a set of finite-difference 
equations can be described by the phase speed of the physical mode. The 
true phase speed of th<^ one-dimensional shallow-water system is (U + C ), 

O 

After thorough analysis, the phase speeds of different approximations in 
space for this system are the same, i.e., (U + C^) among 

the mi:ced, pure, and leap-frog methods. Apparently, the true solution is 
generally retarded by any one of the approximations. In other words, 
the computational dispersion of the mixed method is as good as the leap- 
frog or pure method, 
d. Two-dimens ior-"^il case 

For a two-dimensional shallow-water system, the Coriolis effect may 
be incorporated. The system consists of the folla<ring equations. 


8u 

yt 


+ 


tI + V ^ 


g 


3x 


- fv = 0, 


3v 

3t 


+ 


^ 3x '' 3y 


g 


3h 

5y 


+ fu = 0 


(25. a) 
(25. b) 


3t 

3h 


+ 


U 

3x 


3y 


H (■ 


3u , 3v 

w 


) = 0 


(25. c) 


The mixed method is applied to this system except u and v in the Coriolis 
terms are approximated by u^ ^ and ^ respectively. The two-dimensional 
stability criterion for the difference system is 


At ^ 


B 


(26) 
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31 « 


where 






A = (1-cos cxtx) + 2(^)2 (1-cos Piy) 

+ 2 ^ ^ 1-cos aAx - cos 0Ay + cos(aAx - 3Ay> 

(27) 

■i (f^ + 6^) 

g 

+ 2(^ sin aAx -f ^ sin 3Ay) (f^ + 


B = 2 ~ (1-cos aAx) + (1-cos 3 Ay) 


(28) 


and 


{2 _ sin^ aAx ^ sin^ 3Ay 


Ax" 


Ay 


(29) 


if the solutions have the form 



\ exp l^i (aJlAx + 3mAy)J . 




(30) 


It is clear that the Coriolis effect is dominated by the surface 
gravity waves, if the fluid is deep enough. The basic properties of 
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this two-dimensional approximation is quite similar to the previous 
one-dimensional one. 


3. Linear primitive-equation: A three-dimensional system 

To model atmospheric phencmena in mesoscale, a more complete set 
of coupled primitive equations than a shallow water system is clearly 
needed. Careful analysis of numerical stability for the approximate 
system is necessary to ensure stable calculations. The linear primi- 
tive equations which will be used here are the simplifications of a 
three-dimensional mesoscale model equations (Hsu, (8)). They are 




(31.=’' 




3v . n 9v . r; 3v . Q .. 3 v ^ 


3tt 




(31. b) 
(31. c) 


3x 3y 3z 


(31. d) 


39 

3t 




+ V-||+wS-K-|^=0 


(31. e) 
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a. Stability 

Let 1 denote the finite-difference approximation to 

)w } lU ) iC 

A(t,x,y,z) = A(nAt,ilAx, mAy, IcAz) for the dependent variables u, v, w, 
■n, and e. The finite-difference equations of the mixed method are 


and 


U« - Un 1 U - U T TToj.1 ~ 

Z + n J + V -2 + 00 -iii — Sii 


At 


Ax Ay 

n Vl " Vl _ 


2Ax 


- fv“ - K 


T 


0 , 


Az‘ 

, ij - Vl . r, VVl , 0 , Vl ~ Vl 


At 


Ax 


Ay 


2Ay 


- 0 , 


yn ~ ^Tc 

Az 


Az 

Vi llK.o. 


^Jl+1 *^L-1 ^ ^m+1 ^tn-1 ^ ^k+1 ^k-1 _ q 


0 


2Ax 

n+1 _ gn 


2Ay 


2Az 

0 - e. 


+ u!LlfM + V>— Sw" 


At 


Ax 


Ay 


,, ®k+l ” ^®k ®k-l _ Q 

k 7 u, 

Az 


(32, a) 


(32. b) 


(32. c) 
(32. d) 


(32. e) 


when we assume that 11 > 0, and V > 0. 
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Wave solutions of (32) take the form 

m k “ ^ exp{i(aaAx+BmAy+YkAz)}. (33) 

After substituting (33) into (32) and performing some algebraic manipu- 
lations, an amplification matrix can be obtained, 


where 


G = 


\l)-i4) 

fAt 

-^x 

-fAt 




s 

\}j-i4) 


At At 

= 1 - 1j — (1 - cos aAx) - V — (1 - cos 0Ay) 
Ax Ay 

KAt 

+ 2 — R* (1 - cos yLz ) , 

Az^ 

At At 

(j) = Xy — sin aAx + V ■ — sin pAy,‘ 

Ax Ay 


(34) 


(35) 

(36) 




5„ = 


and 


g At 

Az 

sin aAx 

00 

sin yAz 

Ax , 

g At 

Az 

sin 0Ay 

00 

sin yhz 

4y , 

S At 

Az 

sin aAx 

sin yAz 

Ax , 

S At 

Az 

sin SAy 

sin yAz 

Ay . 


(37) 


(38) 


(39) 


(40) 
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The physical eigenvalue of the amplification matrix is 


X = 4, - i{* + 

The squared absolute value of X is 


(41) 


X| = 1 - BAt + AAt^, 


(42) 


where 


A = 2(^)\l - cos aAx) + 2(^j*(l - cos BAy) 


U V 


+ 2 > 2 ^ {1 - cos aAx - cos 3 Ay + cos (oAx-gAy)} 


+ 2(i sin aAx + X sin gAy)(N^6^ + f^)^ 


+ (N^5^ + f2) 


+ 4(1 “ "^1 - cos aAx) - ~(1 - cos 3Ay)} 


K 

{ — jv (1 “ cos yAz)} 
Az^ 


+ 4 (1 ” cos yAz)^ 

Az^ 


(43) 




( 44 ) 
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and 


6 ^ = 


te‘ 

:2(jL - cos VAz) 


2 . 2 
rSin aAx ^ sin 3Ay ^ 


N‘ 



S d9o 
oqc^z • 


(45) 

(46) 


Then the stability criterion is given in the von Neumann’s sense, 


At £ ^ (47) 

For this three-dimensional primitive -equation system, the mixed 
method is conditionally stable . 

Since |Xi >1. This method becomes absolutely unstable, 

so U or V is not allowed to vanish. 

b. Averaging in the hydrostatic equation 

The averaging procedure appears in (32. c), and is crucial to the 
stability of the approximate system. If 0 in (31. c) has not been 
averaged, we may instead have 


„ Vi " Vi __ s „ 
5S5 ® ‘ 


.n 


(48) 


The only affected term in (42) is (45) and becomes 


2 .2 .2 

p2 _ Az ^ sin aAx , sin 3Ay v 

c, ^ VT ' ; 


. 2 . Ax‘ 

sin yAz 


(49) 


A/ 


2 ^ 

For 2Az i^ve, 6 is unbounded, and |X[ is much greater than 1. Hence 


the unstable situation occurs due to the inproper approximation of the 
hydrostatic equation, 
c. Influences of physical factors 

Generally, the relationship between U (or V) and At in the primi- 
tive-equation systems still remains quite similar to that in the 
shallar^ter system, because the same numerical technique is employed. 
However, two systems describe different physical waves, i.e. surface 
gravity waves for the shallow^-water system and internal gravity waves 
for the primitive-equation system. 

In the present primitive-equation system, three physical factors 
govern the stability bound and they are 

(1) The Coriolis effect (K = N = 0) 

(2) The thermal stratification (§ = K = 0), and 

(3) The vertical diffusion (N = $ = 0) . 

l^en any one of these three factors is retained in the system At will have 
a maximum with respect to U (or V) . They are illustrated in figures 
(4a), (4b), and (4c) for the case (!)-§= 40°, case (2) - N = 0.01 
sec , and case (3) - K = 10*^ cm sec respectively. 

At decreases from its maximum to zero as U (or V) either increases 
to infinite or decreases to zero. These results are different from the 
case which none of the physical factors appears in the system ($ = K = N = 0) . 
In figure (4d), At is inversely proportional to U (or V), and At is 
unbounded as U (or V) goes to zero (i.e. the CFL condition). 

Usually, At increases as Ax (Ay) increases. It can be found in 
figures (4b) and (4d) . Both in figure (4a) and (4c) , At increases as 
Ax (or Ay) decreases under weak U (or V) conditions. It implies that 
under some situations, the shorter the horizonfoil increment, the larger 
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the time increment. This is quite unusual. In applications,' At 
should be determined by exact calculation of (47). 

d. Vertical dependence 

The vertical dependence of the stability bound appears in the 
thermal stratification and vertical diffusion terms in (43). The 
limiting condition of that both the vertical increment and the vertical 
^^avelength approach to zero provides the same result as eliminating the 
thermal stratification (N = 0) and the vertical diffusion (K= 0) 
effects (Fig. 4a). 

If only the vertical diffusion effect is omitted in (43), the max- 
imum time step decreases as the vertical increment and the vertical 
wavelength increase (F~g. 5a) . This means that the shorter the vertical 
increment, the larger the time step. High accuracy in vertical with a 
long time step is obviously allowable. The situation becomes completely 
opposite when only the thermal str*atification effect is not considered 
(Fig. 5b). \Vhile both effects are retained in the primitive equation 
system (Fig. 5c), it seems that the stability bound of a pure thermal 
stratification case (Fig. 5a) is modified by the vertical diffusion 
effect. Furthermore, for longer ^^ves the modification of time increment 
by the vertical diffusion is weaker than that for the shorter waves. 

The reduction of time step is quite significant for short wavelength 
and short vertical increment. 

e . Accuracy 

In general, the accuracy of the mixed method applied to the three- 
dimensional primitive-equation system is quite good. Figure 6 shows 
Che plot of |x|^ vs. At at N = 0.01 sec”^, K = 10^ cm^ sec”^, | = 40°, 

U = V = 10 m sec 10 Ax, 10 Ay and Ax = Ay = 30 Icm. The la-zest value 
of |;\p is 97.63% (i.e. jxj = 98.81%) for the 6Az waves. This very 


50 . 


weak damping may be very helpful to suppress the instability caused by 
non-linear integration, 

4. Conclusion 

A mixed numerical method has been developed for atmospheric models , 
and consists of the forward difference scheme for time tendency terms, 
the upstream scheme for advection terms, and the central scheme for 
other terns in a physical system. For simple advection equation, the 
f onward-upstream method is excessively dissipative , awhile the f onward- 
central one is absolutely unstable . The mixed method is a combination 
of these t\wo methods. Most importantly, this method is conditionally 
stable and highly accurate to the approximate system of either the 
shallow-water equations in one or two dimensions or the primitive equa- 
tions in three dimensions. The dependences of determining the stability 
bounds are not quite obvious. However, the analytic expressions of the 
linear stability criteria are given. At should be easily found under 
typical conditions. Tlie traditional CFL criterion is only a special 
case of the present results, which give opposite criterion to the CFL 
criterion in certain situations, 

The mixed method not only conserves computer resources but also 
programming efforts, because it is explicit and two-tima-level. This 
method has been successfully applied by the author in his mesoscale 
model (Hsu, (8). Stable calculations have been achieved \without any 
artificial spatial smoother or temporal filter . 

To compare the accuracy of the model results among the mixed method 
and other explicit ones for a non-linear atmospheric model, it is nec- 
ess^>ry to simulate some real cases with different methods and analyze 
the differences betsween observational and computational data. This 



effort is being undertaken. Implicit or semi“implicit method (e.g. 
Sun, (9) may be helpful to increase the length of the allowable time 
step for long-term integration. 
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Appendix 


List of Symbols 


C speed of surface gravity wave 

O 

f 2flsin4>, Coriolis parameter 

g gravitational acceleration 

G amplification matrix 

h height perturbation of a free surface 

H constant height of a free surface 

i 

K vertical eddy exchange coefficient 

il, m, k number of increments in x-, y, and z~ directions, respectively 


N 


g 

(— 


de 

0 

dZ 


i 


Brunt-Vdisala frequency 


de 

S constant potential temperature lapse rate 

t time 

u, V, w X-, y-, and z- components of velocity perturbation, respectively 

U, V constant velocities in x~ and y- directions, respectively 

U value of U corresponding to At 

max r o 

X, y, z Cartesian coordinate 

a, 3, Y ^^venumbers in x~, y-, and z- directions, respectively . 

At time increment 

At^y limiting At as the wave number approaching to zero 

^^max maximum of At.^ 


Ax, Ay, Az space increments in x-, y-, and z- directions, respectively 
0 potential temperature perturbation 

0^ constant potential temperature 

X eigenvalue of the amplification matrix 

IT scaled pressure perturbation 

$ latitude 

n Earth rotation rate 
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FIGURE LEGENDS 


Figure 1: 


Figure 2: 


Figure 3: 
Figure 4: 

Figure 5: 


Figure 6: 


Variations of At.., with respect to U for (a) H = Q, (b) H = 1 m, 
(c) H = 10 rn, (d) H = 100 tn, and (e) H = 1000 m. Curve lables 
are Ax in 10 km. 

2 

Damping factor |X| of the mixed method for (a) 2Ax - 4Ax waves 
and (b) 4Ax ~ lOAx waves (Ax = 30 km, U = 10 m-sec and H = 
1000 m) . ^ 

Damping factor \X\ of the pure method for 2Ax - lOAx waves 

( X = 30 km, U - 10 mrsec ^ , and H = 1000 m) . 

Variations of At with respect to either U or V for (a) K = 0, 

N s 0, 40°, (b) K = 0, N = 0.01 sec"^, $ = 0, (c) K = 10^ 

cm -sec N = 0, <I> = 0, and (d) K = 0, N = 0, <I> =0 for lOAx- 

10Ay-5Az wave (Az = 1 km). Curve lables are Ax=Ay in km. 

Variations of At with different vertical (Az) waves for (a) 

K = 0, N = 0.01 sec”^, (b) K = 10^ cm^-sec~^, N = 0, and (c) 

•3 9 , _1 -il , 

K = 10 cm -sec , N = 0.01 sec . Horizontal constants are 

^ 

10Ax=10Ay waves with Ax=Ay= 30 km, $ ■= 40° , and U=V= 10 ra-sec . 

Curve lables are Az in l<m. 

2 

Damping factor [X]- of the mixed method in the three-dimensional 
mesoscale model for different vertical (Az) waves. Constants 
are 10Ax=10Ay waves with Ax=Ay = 30 km, U = V = 10 cm-sec 
$ “ 40°, K = 10^ cm^-sec \ and N = 0.01 sec 
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